Our project will try to take advantage of machine learning methods to try to predict, with ~100 ms in advance, the onset of a “seizure". The recordings used in the development phase were from spiking activity of cortical neurons on microelectrode arrays chips with 256 electrodes. These recordings exhibit seizure-like activity, where a large fraction of the population synchronizes.
Epileptic seizure is a period of symptoms due to abnormally excessive or synchronous neuronal activity in the brain.
First some modules need to be imported:
# These are the imports of the McsData module
import sys
sys.path.append('D:\\Programming\\McsDataManagement\\McsPyDataTools\\McsPyDataTools')
import McsPy.McsData
import McsPy.McsCMOS
# numpy is numpy ...
import numpy as np
# bokeh adds more interactivity to the plots within notebooks. Adds toolbar at the top-right corner of the plot.
# Allows zooming, panning and saving of the plot
import bokeh.io
from bokeh.io import show
from bokeh.layouts import column
from bokeh.models import ColumnDataSource
from bokeh.plotting import figure
from bokeh.models.widgets import Slider
bokeh.io.output_notebook()
# tensorflow & other machine learning utilities
import tensorflow as tf
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import normalize
# utilities
import time
Next we need to access the rawdata by initializing an instance of the RawData class from the McsData module by handing over the path to the file. The filepath points to the folder TestData within the folder, where this notebook resides.
To check if we got access to the file we can look at its contents by printing the info that got extracted when the RawData object was initialized.
channel_raw_data = McsPy.McsData.RawData('./TestData/cortical_same_chip/2019-02-11T11-42-58Chip_3_DIV17.h5')
print(channel_raw_data.comment)
print(channel_raw_data.date)
print(channel_raw_data.clr_date)
print(channel_raw_data.date_in_clr_ticks)
print(channel_raw_data.file_guid)
print(channel_raw_data.mea_name)
print(channel_raw_data.mea_sn)
print(channel_raw_data.mea_layout)
print(channel_raw_data.program_name)
print(channel_raw_data.program_version)
The data of the stream can be found under .channel_data.
analog_stream_0 = channel_raw_data.recordings[0].analog_streams[0]
analog_stream_0_data = analog_stream_0.channel_data
#analog_stream_0_data = analog_stream_0_data[:,200000:300000]
electrodes = analog_stream_0.channel_data.shape[0]
print(analog_stream_0_data)
def modify_doc(doc):
electrode_index = 0
source = ColumnDataSource(data=dict(x=range(analog_stream_0_data.shape[1]), y=analog_stream_0_data[electrode_index]))
bfig = figure(plot_width=900, plot_height=400, title='Voltage Activity (volt) / Time (microsecond)')
bfig.circle('x', 'y', source=source)
bfig.xaxis.axis_label = 'Microsecond'
bfig.yaxis.axis_label = 'Voltage'
bfig.ygrid.minor_grid_line_color = 'navy'
bfig.ygrid.minor_grid_line_alpha = 0.1
def callback(attr, old, new):
source.data = ColumnDataSource(data=dict(x=range(analog_stream_0_data.shape[1]), y=analog_stream_0_data[new-1])).data
slider = Slider(start=1, end=electrodes, value=1, step=1, title="electrode")
slider.on_change('value', callback)
doc.add_root(column(slider, bfig))
show(modify_doc) #show(modify_doc, notebook_url="localhost:8890")
The detection of spikes will follow the formula:
spike = abs(voltage) > abs(std_electrode_voltage*factor)
You can define the factor below (using 5 as default)
def aux(value, std, factor):
return False if abs(value) < abs(std*factor) else True
def threshold(electrode_values, factor=5):
std = np.std(electrode_values)
return np.array([aux(value, std, factor) for value in electrode_values])
start = time.time()
np_analog_stream_0_data_spikes = np.array([threshold(electrode_values, 5) for electrode_values in analog_stream_0_data])
end = time.time()
print("Elapsed time: " + str(end-start))
print(np_analog_stream_0_data_spikes)
def modify_doc(doc):
electrode_index = 0
colors_spikes = ['green' if spike else 'blue' for spike in np_analog_stream_0_data_spikes[electrode_index]]
source = ColumnDataSource(data=dict(x=range(analog_stream_0_data.shape[1]), y=analog_stream_0_data[electrode_index], color=colors_spikes))
bfig = figure(plot_width=900, plot_height=400, title='Voltage Activity (volt) / Time (microsecond)')
bfig.circle(x='x', y='y', source=source, color='color')
bfig.xaxis.axis_label = 'Microsecond'
bfig.yaxis.axis_label = 'Voltage'
bfig.ygrid.minor_grid_line_color = 'navy'
bfig.ygrid.minor_grid_line_alpha = 0.1
def callback(attr, old, new):
colors_spikes = ['green' if spike else 'blue' for spike in np_analog_stream_0_data_spikes[new-1]]
source.data = ColumnDataSource(data=dict(x=range(analog_stream_0_data.shape[1]), y=analog_stream_0_data[new-1], color=colors_spikes)).data
slider = Slider(start=1, end=electrodes, value=1, step=1, title="electrode")
slider.on_change('value', callback)
doc.add_root(column(slider, bfig))
show(modify_doc) #show(modify_doc, notebook_url="localhost:8890")
The detection of seizures will follow the formula:
seizures = #electrodes spiking >= factor * electrodes
You can define the factor below (using 2.5% as default)
spike_times=[]
electrode_spike=[]
seizure_dict={}
for electrode, line in enumerate(np_analog_stream_0_data_spikes):
for time, val in enumerate(line):
if val:
spike_times.append(time)
electrode_spike.append(electrode)
set_electrodes = seizure_dict.get(time)
if set_electrodes is not None:
set_electrodes.add(electrode)
seizure_dict[time] = set_electrodes
else:
seizure_dict[time] = {electrode}
factor=0.05
threshold = electrodes * factor
seizure_detection = np.zeros(np_analog_stream_0_data_spikes.shape[1], dtype=bool)
for time in seizure_dict:
spikes = len(seizure_dict[time])
if spikes > threshold:
seizure_detection[time]=True
colors = ['red' if seizure_detection[t] else 'blue' for t in spike_times]
plot_size_and_tools = {'plot_height': 400, 'plot_width': 900,
'x_range': [0,np_analog_stream_0_data_spikes.shape[1]], 'y_range': [0, np_analog_stream_0_data_spikes.shape[0]]}
raster = figure(title="Raster Plot", **plot_size_and_tools)
raster.circle(x=spike_times, y=electrode_spike, color=colors)
show(raster) #show(raster, notebook_url="localhost:8890")